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Abstract  »— J — — 

^  The  world/of  mathematical  programming  has  seen  a  remarkable  surge  of 
activity  following  publication  of  Karmarkar’s  projective  alg^itluQjfi_May  _^984. 

A  review  of  me  ensuing  three  yeairs  has  already  appeared:  (Montiia  "f34})^One 
year  later,  we  review  some  of  the  main  methods  eind  surrounding  events,  and 
focus  on  references  that  contain  computational  results.  \ 

C  ^ -  '  'X  .  ,  . 

Keywords:  Linear  programming,  interior-point  methods,  barrier-function 
methods,  Newton’s  method,  Karmarkar’s  projective  method^  ^  , 

1.  Introduction 

It  is  still  only  four  years  since  Narendra  Karmarkar  of  AT&T  Bell  Laboratories 
presented  the  mathematical  features  of  an  apparently  new  method  for  solving  linear 
programming  (LP)  problems  [23].  The  problem  was  assumed  to  be  in  the  form 
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minimize 


subject  to  Ax  =  0,  e^x  =  1,  a;  >  0, 


where  A  is  m  x  n,  m  <  n  and  e  is  a  vector  of  I’s.  It  was  also  assumed  that  x  =  e 
is  a  feasible  solution,  and  that  the  optimal  objective  value  is  zero.  These  and  other 
restrictions  have  since  been  lifted  (though  not  without  practical  difficulty). 

*The  material  contained  in  this  report  is  based  upon  research  supported  by  the  Air  Force  Office 
of  Scientific  Research  Grant  87-01962;  the  U.S.  Department  of  Energy  Grant  DE-FG03-87ER25030; 
National  Science  Foundation  Grants  CCR-8413211  and  ECS-8715153;  and  the  Office  of  Na.al  Re 
search  Contract  N00ni4-87  K-(il42. 
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A  vital  part  of  Karmarkar’s  analysis  involved  showing  that  the  so-called  potential 
function  c^xj  could  be  reduced  by  a  fixed  amount  every  iteration.  It  followed 
that  even  in  the  worst  case,  the  number  of  iterations  required  to  reach  a  solution 
would  be  proportional  n.  The  work  per  iteration  wa.s  shown  to  be  of  order 
Hence,  the  total  time  required  was  bounded  by  a  polynomial  in  n. 

Broadly  speaking,  the  bound  on  computation  time  for  the  projective  method  is 
of  order  (where  T  is  a  measure  of  the  storage  needed  for  ^he  input  data).  In 

more  recent  methods  the  bound  has  been  reduced  to  0{n^L)  [34,37],  whereas  for 
the  classical  simplex  method  [4]  the  total  time  on  certain  contrived  problems  is  of 
order  n^2”. 

1.1,  A  practical  algorithm? 

To  operations  researchers  and  others,  it  seemed  immediately  evident  that  the  amount 
of  work  per  iteration  for  the  new  method  would  be  much  greater  than  for  the  simplex 
method,  at  least  for  “normal”  LP  problems.  The  crucial  part  of  the  computation 
involves  a  Cholesky  factorization  AX^AJ  =■  R^R  (with  X  diagonal,  R  triangular), 
and  this  is  normally  much  more  expensive  in  time  and  storage  than  the  sparse  basis 
factorization  B  —  LU  required  by  the  simplex  method. 

Nevertheless,  it  seemed  that  the  reduced  number  of  iterations  might  sometimes 
compensate,  and  for  very  specially  structured  problems  it  was  reasonable  to  suppose 
that  the  Cholesky  factors  could  be  obtained  quickly  enough  to  give  an  efficient 
algorithm. 

1.2.  Not  just  practical 

Controversy  soon  enveloped  the  new  algorithm.  In  a  series  of  conference  appearances 
that  are  continuing  to  this  day.  Karmarkar  has  reported  results  that  are  uniformly 
superior  to  those  obtained  by  the  simplex  method.  The  speed-up  factors  claimed 
vary  from  4  up  to  84,  104,  190,  720  and  even  higher  [24]. 

The  factor  of  4  was  obtained  using  part  of  a  set  of  13  publicly  available  problems, 
supplied  on  request  by  Stanford  University  in  July  1984.  This  test  set  was  later  dis¬ 
missed  as  “not  representative  of  modern  real-world  problems  in  size  or  complexity”, 
even  though  several  problems  were  omitted  in  the  results  reported,  including  the 
problem  that  was  by  far  the  largest!  (The  PILOT  model  had  about  1500  rows,  3700 
columns  and  43000  nonzeros — certainly  only  medium-scale  by  conventional  stan¬ 
dards,  yet  “large”  in  the  sense  that  a  cold-start  solution  with  the  simplex  method 
takes  over  20  hours  on  a  DEC  VAXstation  II.') 

Although  much  larger  problems  are  indeed  of  interest,  credibility  at  the  time 
depended  strongly  on  results  for  lamiiiai  ;noucls  Omission  rvithout  comment  of  the 
known  most-difficult  test  cases  inevitably  generated  scepticism. 

Further  controversy  arose  from  the  fact  that  the  more  spectacular  results  were 
obi,<ti.t..u  fror.i  piuLLii.o  thci,'.  proprietary,  like  the  implementation  itself.  Nei¬ 
ther  the  problems  nor  the  program  could  be  seen  by  the  outside  world,  and  apart 

’This  is  a  typical  run-time  for  MINOS  5.3  (May  1988),  a  portable  Fortran  code.  A  commercial 
.Mathematical  Programming  System  would  lake  about  30  minutes  on  an  IBM  3090. 
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from  problem  dimensions,  the  only  results  given  were  iteration  counts  and  CPU 
time.  No  plot  of  the  individual  iterations,  no  stopping  criteria,  no  details  about  the 
precision  obtained  or  the  reliability  of  the  code. 

1.3.  A  flurry  of  effort 

In  spite  of  the  unsatisfactory  circumstances,  much  feverish  activity  has  been  gener¬ 
ated  within  the  scientific  community,  and  numerous  researchers  have  launched  into 
projects  they  did  not  expect  to  be  pursuing. 

This  is  not  entirely  due  to  a  windfall  of  government  funding!  There  is  also  the 
intrigue  and  the  challenge.^ 

1.4.  An  objective  view 

The  connections  with  more  traditional  areas  of  nonlinear  programming  and  nu¬ 
merical  linear  algebra,  along  with  much  analysis  of  path-following  methods  (e.g., 
[29,30,35])  have  cast  sufficient  light  on  the  scene  for  us  to  believe  that  good  perfor¬ 
mance  is  indeed  possible  on  a  significant  proportion  of  real-world  problems. 

In  terms  of  robustness  the  verdict  is  still  out,  since  present  implementations 
(within  our  experience)  are  highly  sensitive  to  slight  changes  in  strategy.  However, 
it  took  four  decades  to  make  the  simplex  method  (virtually)  100%  reliable.  Much 
has  been  learned  in  that  time  about  linear  and  nonlinear  optimization.  While  the 
methods  inspired  by  Karmarkar  are  certainly  difficult  to  make  “bullet-proof’,  there 
is  every  hope  that  useful  and  acceptably  reliable  implementations  will  be  developed 
within  the  next  few  years. 

In  the  following  sections  we  review  the  main  mathematical  variants  that  have 
been  proposed  so  far,  giving  known  computational  results  where  possible. 

2.  Primal  and  Dual  LPs 

The  special  form  (1)  allowed  certain  geometrical  arguments,  and  many  research 
papers  have  been  written  around  it.  However,  operations  researchers  have  long 
been  generating  linear  programs  in  the  more  practical  form 

Primal:  minimize  c^x  subject  to  Ax  =  b,  x  >  0,  (2) 

with  no  assumptions  about  the  optimal  objective  value.  Variants  of  the  Karmarkar 
approach  began  to  appear  for  solving  (2)  directly  (without  converting  it  to  the 
special  form);  e.g.,  [12,18].  It  was  soon  found  that  dual  variables  could  be  obtained 
as  by-products.  Algorithms  for  the  dual  problem 

Dual:  maximize  b^n  subject  to  <  c,  (3) 

^The  simplex  method  has  been  extremely  successful  so  far,  and  computers  are  growing  ever 
faster,  yet  it  is  never  enough.  Factors  of  100  or  more  are  tantalizing  and  compelling,  even  in  the 
face  of  scepticism. 
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also  began  to  appear.  We  discuss  the  deveiopments  for  both  problems  in  roughly 
chronological  order.  The  most  general  case  with  upper  and  lower  bounds  on  x  is 
covered  in  Section  6. 
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3.  Newton’s  Method  and  Barrier  Functions 

Given  a  vector  ar,  let  A"  be  a  diagonal  matrix  whose  y-th  diagonal  is  xj.  It  will  be 
assumed  that  xj  >  0  for  y  =  1  to  n.  Similarly  for  quantities  z  and  Z. 

The  barrier  function  approach  deals  with  inequality  constraints  such  as  x  > 
0  by  adding  a  judicious  function  to  the  true  objective.  A  characteristic  is  that 
along  any  direction  of  descent  from  an  interior  point,  there  exists  a  strictly  interior 
minimizer.  Usually  the  transformed  objective  function  is  infinite  along  the  boundary 
of  the  feasible  region.  Common  barrier  functions  are  the  logarithmic  and  reciprocal 
functions. 

3.1.  A  primal  barrier  method 

For  the  primal  LP  we  consider  subproblems  of  the  form 

n 

minimize  F^{x)  =  c^x  -  fi^]nxj  subject  to  Ax  =  b,  (4) 

*  i=i 

where  /x  (^  >  0)  is  a  specified  parameter  that  will  be  set  to  decreasing  values.  As 
p  0,  the  solution  of  (4)  approaches  that  of  (2). 

One  of  the  first  computational  studies  outside  AT&T  was  undertaken  by  Gill  et 
al.  [18],  who  recognized  the  connection  between  Karmarkar’s  method  and  Newton’s 
method  applied  to  the  barrier  subproblems.  Newton’s  method  obtains  a  search 
direction  p  by  solving  a  system  of  the  form 


where  g  and  H  are  the  first  and  second  derivatives  of  F^i{x). 

Given  a  feasible  interior  point  *,  an  estimate  of  the  dual  variables  tt  and  the 
corresponding  “reduced  costs”  z  =  c  —  A^n,  the  main  steps  prove  to  be  as  follows: 

1.  Define  r  =  z  -  pX~^e. 


2.  Solve  the  system 

AX\\^q  =  AX^r,  or  equivalently,  min  ||A(r  -  A^7)||2.  (5) 

3.  Update  x-  ♦—  tt  f  q. 

4.  Define  z  =  c  -  A^ir,  r  =  z  -  pX~^€  and  p  =  -(l//f)A^r. 

5.  Find  a  steplength  a  at  which  the  barrier  function  +  ap)  is  suitably  less 
than  Ff,(x). 
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6.  Update  a;  <—  i  +  ap. 

The  algebra  is  evidently  not  much  more  complicated  than  in  the  simplex  method. 
However,  Step  2  is  crucial.  The  matrix  AX'^A^  is  large,  it  could  be  quite  dense 
compared  to  A,  and  in  general  it  is  very  nearly  singular.  A  combination  of  direct 
and  iterative  methods  may  be  applied,  but  the  system  must  be  solved  accurately  in 
order  to  retain  the  condition  Ax  =  6. 

If  the  barrier  parameter  p  is  held  constant,  the  vector  r  will  eventually  become 
small  as  Steps  1-6  are  repeated.  The  iterations  then  start  again  with  a  smaller  p. 
In  practice,  we  typically  give  p  the  values  0.1,  0.01,  . . . ,  10“®  (say)  over  a  total  of 
20-50  iterations. 


3.2.  Relationship  to  Karmarkar’s  projective  method 

The  best  choice  of  p  remains  open  to  question.  In  [18]  it  was  shown  that  if  the 
Newton  barrier  method  is  applied  to  the  special  LP  problem  (1),  and  if  p  and  a 
are  chosen  in  a  certain  (easily  computed)  way,  then  the  sequence  of  points  x  are 
identical  to  those  generated  by  Karmarkar’s  original  method.  In  other  words,  the 
“new”  method  for  linear  programming  may  be  viewed  as  a  special  case  of  Newton’s 
method  for  linearly  constrained  optimization. 

This  is  not  to  minimize  Karmarkar’s  contribution.  For  example,  the  “special 
case”  argument  was  deficient  without  a  proof  that  the  special  choice  of  p  would  be 
positive. 

One  practical  contribution  has  been  Karmarkar’s  unrelenting  advocacy  of  the 
approach,  even  though  it  requires  solution  of  systems  involving  AX^A^.  Years  earlier, 
prior  to  modern  sparse-matrix  techniques  (e.g.,  [14]),  researchers  had  rejected  similar 
approaches  as  being  unacceptably  slow  on  large  problems. 

In  addition,  Karmarkar’s  particular  choice  of  fi  and  a  enabled  him  to  obtain 
a  polynomial  bound  on  the  iterations  and  work  required.  (No  such  bound  was 
previously  known  for  barrier  methods.)  Recently,  Gonzaga  [22]  has  shown  that  when 
Karmarkar’s  method  is  generalized  to  treat  problem  (2)  it  again  corresponds  to  the 
barrier  method  with  a  certain  positive  choice  of  p.  Hence,  the  projective  method 
and  the  primal  barrier  method  are  now  known  to  be  polynomial-time  algorithms 
even  on  the  more  useful  formulation  (2). 

Note  that  the  Newton-barrier  approach  has  practical  and  theoretical  benefits. 
In  contrast  to  Karmarkar’s  original  algorithm,  the  extension  to  other  forms  of  LP 
and  to  nonlinear  programs  is  trivial.  The  Newton-barrier  approach  is  also  the  basis 
of  the  analysis  in  [34,37]. 

3.3.  Results 

In  [18],  results  were  given  for  nine  of  the  problems  sent  to  Karmarkar  (see  Sec¬ 
tion  1.2);  the  remainder  had  upper  bounds  and  could  not  be  run.  In  eight  of  the 
nine  cases,  the  primal  barrier  method  was  found  to  be  comparable  in  speed  to  the 
simplex  method,  and  on  two  additional  (larger)  models  a  speed-up  factor  of  3  was 
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obtained.^'** 

Unfortunately,  a  further  series  of  highly  degenerate  models  showed  an  unfavor¬ 
able  trend  with  increasing  problem  size.  Nevertheless,  the  conclusions  in  [18]  struck 
an  optimistic  note  for  interior-point  methods,  particularly  for  problems  with  special 
structure. 

More  recently,  a  new  primal  barrier  implementation  has  been  developed  at  Stan¬ 
ford  [16].  Of  particular  interest  is  the  solution  time  for  PILOT:  about  9  hours  on 
a  VAXstation  II.  This  is  a  speed-up  of  2.3  on  a  real-world  model  that  is  unques¬ 
tionably  non-trivial  for  the  simplex  method.  The  periodic  structure  revealed  in  [26] 
may  be  a  contributing  factor,  but  in  any  event,  this  represents  a  bright  note  for 
the  interior-point  approach  within  the  scope  of  current  repeatable  computational 
results. 

3.4.  A  dual  barrier  method 

The  barrier  subproblem  corresponding  to  the  dual  LP  (3)  is 

n 

maximize  b^ir  -b  //  ~  “J^)» 

J=i 

which  is  a  purely  unconsiramed problem.  For  certain  x  and  X  [19],  Newton’s  method 
leads  to  the  equations 

AX^A^q  =  n{Ax  -  6),  tt  ♦—  tt  +  aq,  (6) 

and  we  can  capitalize  on  the  fact  that  the  system  for  q  need  not  be  solved  exactly  [5]. 
This  opens  to  door  to  obtaining  “cheap”  approximate  factors  of  the  matrix  AX^A^^ 
for  use  as  preconditioners  with  the  conjugate-gradient  method. 

3.5.  Results 

For  reasons  described  by  Gill  et  al.  [19],  a  sparse  factorization  XA^  =  LU  should 
provide  a  useful  preconditioner.  An  experimental  implementation  has  so  far  proved 
to  be  less  efficient  than  hoped.  However,  we  anticipate  that  the  advantages  of 
approximate  factors  will  ultimately  come  to  the  fore  (see  Section  7.3). 

4.  Affine  Scaling 

Independently  of  the  barrier-function  development  described  above,  groups  were  at 
work  on  algorithms  for  the  primal  problem  (2)  and  the  dual  problem  (3).  These 
were  the  so-called  affine  variants  of  Karmarkar’s  method,  which  seemed  at  least  to 
be  simpler  than  the  projective  method. 

^The  results  reported  here  and  in  later  sections  were  obtained  using  “cold  starts”.  The  ability  to 
utilize  “good”  starting  information  remains  virtually  unique  to  the  simplex  method.  The  entropy- 
based  method  of  Eriksson  and  the  shifted  barrier  method  also  have  this  advantage;  see  Section  8. 

^Comparisons  with  the  simplex  method  were  made  using  various  versions  of  MINOS  [36]  without 
scaling  or  partial  pricing.  For  future  experiments  we  would  in  general  recommend  specifying  SCALE 
OPTIOII  2  and  PARTIAL  PRICE  10;  see  Lustig  [26]. 
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4.1.  A  primal  affine  method 

At  Bell  Labs,  Vanderbei  et  al.  [40]  developed  an  algorithm  that  proved  to  be  the 
limiting  case  of  the  primal  barrier  method  as  /x  — >  0  [18].  It  was  subsequently  found 
that  the  method  was  first  proposed  by  Dikin  in  1967  [6]. 

Search  directions  are  generated  by  the  same  system  (5).  Any  hint  of  quadratic 
convergence  (associated  with  Newton’s  method)  is  lost  by  taking  /x  =  0,  and  to  date 
there  has  been  no  proof  of  polynomial-time  complexity.  Nevertheless,  respectable 
computational  results  were  obtained  from  a  production  code  developed  at  Bell  Labs 
(Chen  [3]). 

We  believe  that  a  primal  affine  algorithm  is  embodied  within  a  combined  soft- 
waxe/hardware  system  that  is  currently  being  marketed  by  AT&T:  the  KORBX^'^ 
Linear  Programming  System,  which  includes  an  Alliant  multi-processor  computer. 

4.2.  Dual  affine  methods 

Under  the  guidance  of  Karmarkar,  researchers  at  Berkeley  [2]  developed  an  analo¬ 
gous  dual  algorithm  that  again  proved  to  be  the  limiting  case  of  the  (dual)  barrier 
method  as  /x  —►  0.  Several  advanced  implementation  techniques  were  employed  [1], 
and  promising  results  were  obtained  (see  next  section). 

A  similar  implementation  was  developed  in  1986  at  Bell  Communications  Re¬ 
search  by  Monma  and  Morton  [32]. 

4.3.  Results 

There  is  now  a  collection  of  over  50  test  problems  available  publicly  via  netlib  (see 
Gay  [10],  Lustig  [26]).  About  30  of  these  problems  were  used  to  test  the  dual 
affine  implementations  just  mentioned.  Both  codes  achieved  average  speed-ups  of  3 
relative  to  MINOS  (versions  4.0  and  5.0  repectively). 

As  before,  some  of  the  more  interesting  problems  were  not  tested  because  the 
implementations  could  not  handle  upper  bounds.  We  remark  that  the  algorith¬ 
mic  parameters  needed  for  satisfactory  performance  on  all  of  the  problems  would 
probably  give  poorer  average  performance  on  the  30-problem  subset. 

Until  recently  it  was  believed  that  a  large  speed  advantage  wa.s  arising  in  the 
Berkeley  implementation  [2]  from  the  use  of  in-line  code  (and  large  amounts  of 
memory)  during  the  Cholesky  factorization  [1].  However,  Gay  [13]  has  now  shown 
that  such  an  advantage  need  not  exist. 


I 


5.  Primal-Dual  Methods 

Barrier  functions  have  recently  been  taken  up  by  several  authors  in  a  primal-dual 
context  (e.g.,  see  Monteiro  and  Adler  [33],  Lustig  [27]).  Briefly,  these  seek  to  solve 
the  primal  and  dual  barrier  subproblems  (4),  (6)  simultaneously.  The  main  equation 
to  be  solved  has  the  same  form  (5)  as  for  the  primal  barrier  method,  except  that 
the  diagonal  matrix  is  replaced  by  XZ~^. 
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The  equations  Ax  =  6,  A^ir  +  z  =  c,  x>0,  z>0  are  satisfied  throughout,  and 
the  iterations  work  towards  satisfying  the  complementarity  condition  XZ  =  n I  as 

/i  — ►  0. 

5.1.  Results 

At  first  sight,  the  primal-dual  approach  would  seem  to  incur  the  disadvantages  of 
both  primal  and  dual  algorithms,  at  least  in  terms  of  obtaining  initial  interior  points. 

Nevertheless,  McShane  et  al.  [28]  have  given  computational  results  for  a  primal- 
dual  implementation,  as  well  as  for  the  dual  affine  code  of  [32]  (evidently  further 
refined).  On  the  30-problem  subset  mentioned  above,  both  codes  showed  an  average 
speed-up  of  about  4  relative  to  MINOS  5.0. 

6.  A  Single-Phase  Dual  Barrier  Method 

In  order  to  run  all  of  the  test  problems  in  the  netlib  collection  it  is  necessary  to 
develop  an  algorithm  that  treats  the  primal  LP  problem  in  its  most  general  form; 

Primal:  minimize  c^x  subject  to  Ax  =  b,  I  <  x  <  u.  (7) 

X 

The  corresponding  dual  problem  may  be  written 

Dual:  maximize  b^z  -  vTy  +  l^z 

jr.y.j  (8) 

subject  to  A^z  —  y  -t-  2  =  c,  y,  ^  >  0. 

An  unexpected  advantage  arising  from  this  formulation  is  that  the  dual  constraints 
A^;r  —  y  -b  z  =  c  are  easily  satisfied  by  a  suitable  choice  of  the  slacks  y  and  z. 
Further  details  are  given  i.i  Gill  et  al.  [17]. 

6.1.  Results 

An  implementation  to  solve  (8)  is  currently  under  development  at  Stanford.  In 
preliminary  tests,  speed-up  factors  in  the  range  1  to  5  relative  to  MINOS  5.2  have 
been  obtained  for  most  of  the  first  53  problems  in  netlib. 

These  results  are  similar  to  those  obtained  for  the  primal  barrier  implementa¬ 
tion  [16]  mentioned  in  Section  3.3,  which  is  designed  for  problem  (7).  Both  codes 
use  essentially  exact  Cholesky  factors  of  AX^A^,  with  iterative  refinement  and  a 
partitioning  (Schur-complement)  scheme  to  handle  dense  columns  of  A.  It  is  hoped 
that  the  dual  implementation  will  show  an  advantage  with  inexact  factorization  and 
preconditioning  (Section  7.3). 

7.  Computational  Matters 

The  issues  of  starting  points,  rank-deficiency,  termination  criteria,  restarts,  etc.  are 
too  lengthy  to  discuss  here.  We  choose  just  four  topics. 
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7.  Computational  Matters 


7.1.  Norma!  equations  versus  least  squares 

In  all  of  the  interior- point  algorithms,  the  search  direction  is  obtained  from  a  system 
of  the  form 

AX^A'^q  =  V,  (9) 

where  X  is  a  diagonal  matrix  and  v  depends  on  the  algorithm.  If  the  right-hand 
side  happens  to  be  of  the  form  v  =  AX^r^  this  system  is  a  set  of  “normal  equations” 
equivalent  to  the  linear  least-squares  problem 


min||X(r  -  A'^q)\\2. 


(10) 


Least-squares  prob'  )ms  can  be  solved  more  reliably  if  treated  as  such.  For  example, 
conjugate-gradient  methods  generally  require  more  iterations  to  solve  (9)  than  they 
do  to  solve  (10),  particularly  when  the  matrix  X^d^is  ill-conditioned  (as  it  invariably 
is  in  this  context). 

We  note  that  not  all  interior-point  methods  permit  the  least-squares  formulation. 
In  this  respect,  the  advantage  lies  with  the  primal  and  primal-dual  variants,  and 
with  the  single-phase  dual  algorithm. 

To  date,  all  implementations  except  [18,19]  have  used  the  (less  reliable)  normal- 
equations  approach. 

7.2.  Boundedness  of  the  projections 

A  clue  to  the  survival  of  interior-point  implementations  in  the  face  of  extreme  iU- 
conditioning  of  XA^has  recently  been  provided  by  Stewart  [39].  In  the  notation  of 
equation  (10),  Stewart  shows  that  as  long  as  A"  is  diagonal  with  nonzero  elements 
and  A  has  full  row  rank,  the  weighted  pseudo-inverse  =  (AX^A^)~^ AX^  and 
the  oblique  projection  A^A^  are  both  bounded  in  norm,  independently  of  X. 

Briefly,  this  means  that  q  and  A^q  in  (10)  will  not  be  arbitrarily  large,  even  if 
X  causes  AA^A^  to  be  almost  singular. 

In  practice,  of  course,  Cholesky  factorization  may  fail  if  AA^A^is  almost  singu¬ 
lar,  and  other  factorizations  may  be  needed;  see  Section  7.4. 

7.3.  Preconditioning 

In  the  linear  programming  context,  virtually  all  implementors  have  used  exact 
Cholesky  factors  of  AA^A^  (excluding  perhaps  a  few  dense  columns  of  A).  The 
main  reason  is  that  the  sparsity  pattern  of  the  normal-equations  matrix  does  not 
change  as  A  changes;  hence  a  single  “analyze”  can  be  performed  on  the  sparsity 
pattern  of  AA^  to  obtain  an  ordering  of  the  rows  of  A  that  preserves  the  sparsity 
of  the  Cholesky  factor.  The  same  ordering  is  used  for  all  subsequent  factorizations 
AA2A'^=  R'^R. 

An  exception  is  Karmarkar  himself.  In  talks  such  as  [24]  he  advocates  forming 
the  normal-equations  matrix  only  approximately: 

M  =  AA'^A^-t-  Eu 
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and  then  factoring  M  approximately: 

N{M  Aal)N'^  =  /  +  £2, 

where  o  is  an  undocumented  spectral  shift,  and  E\  and  E2  are  “small”  matrices 
representing  the  error  involved.  The  matrix  N  is  then  used  as  a  prcconditionor 
when  the  conjugate-gradient  method  is  applied  to  (9).  (It  could  also  be  used  with 
the  lea.st- squares  formulation  (10).) 

A  further  device  is  to  retain  the  same  preconditioner  N  for  seveicJ  iterations 
of  the  main  algorithm  (perhaps  modifying  N  slightly  but  cheaply  at  each  iteration 
[23,38]). 

Are  these  the  keys  to  Karmarkar’s  spectacular  run-times?  Perhaps  so,  at  least 
for  some  problems.  Preconditioning  of  this  kind  is  widely  used  in  other  areas  and 
has  been  studied  in  the  LP  context  by  Gay  (11).  Further  study  is  well  justified,  as 
it  is  for  LU  preconditioning  (Section  3.5).  We  note  that  as  conveigence  occurs,  the 
preconditioner  must  become  more  exact.  For  general  problems,  this  can  be  achieved 
easily  with  LU  factors,  but  not  necessarily  with  Cholesky  factors. 

7.4.  Alternative  factorizations 

While  Cholesky  factors  of  AA'^A^have  been  preferred  to  date  on  efficiency  grounds, 
we  mention  that  greater  reliability  is  likely  with  orthogonal  and  symmetric-indefinite 
factorizations  of  the  form 

and  i)=V^DV. 

For  the  case  where  A  is  sparse,  considerable  progress  has  been  made  recently  by 
George,  Liu  and  Ng  [15]  and  Duff  [7].  Dense  columns  in  A  remain  a  complication 
for  some  problems.  Here  the  dual  algorithms  have  an  advantage  in  not  requiring  a 
(dense)  artificial  column  for  Phase  1. 

8.  Promising  Approaches 

Before  closing,  we  review  some  methods  that  have  yet  to  be  proven  computationally 
but  appear  to  have  significant  computational  advantages. 

8.1.  Eriksson’s  algorithm 

The  primal-dual  approach  was  discu.ssed  in  Section  5.  A  significantly  dijjeixnl 
primal-dual  algorithm  was  given  in  1981  by  Eriksson  [8]  and  further  developed  in  [9]. 
A  primal  subproblem  is  used,  based  on  an  entropy  function  and  a  current  estimate 
of  X  (i*): 

n 

minimize  c^x  —  p  ]n(xj/x'^)  -  {Xj  -  a'j  )| 

2=1 

subject  to  Ax  =  b. 


Primal: 


(11) 
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It  is  assumed  that  x*  >  0,  but  not  that  ylx^  =  b.  Instead,  the  dual  of  (11)  is  treated 
as  an  unconstrained  problem  in  tt: 

n 

Dual:  maximize  6^7r  —  ^  ^  Xje~*j/^,  (12) 

i=i 

where  z  =  c  —  A^ir.  An  inexact  Newton  method  is  applied  to  (12),  with  the  central 
system  of  equations  taking  the  form 

AXA^q-b-Ax'^.  (13) 

This  algorithm  has  many  intriguing  properties,  and  we  believe  it  to  be  of  great 
promise.  For  example,  the  matrix  AX will  in  general  be  better-conditioned  than 
the  usual  AX^A^.  Competitive  computational  results  await  implementation  of  a 
sparse  preconditioner  for  (13),  using  techniques  that  have  been  applied  to  barrier 
algorithms  elsewhere. 

8.2.  The  Box  Method 

One  interpretation  of  the  affine  invariant  algorithms  is  that  they  optimize  the  objec¬ 
tive  function  subject  to  feasibility  and  a  quadratic  constraint.  The  ellipsoid  defined 
by  the  quadratic  constraint  may  be  replaced  by  a  “box”. 

Such  an  approach  has  been  described  by  Zikan  and  Cottle  [41,42],  and  promising 
computational  results  have  been  obtained  for  the  special  case  of  network  problems. 
For  more  general  linear  programs,  an  appealing  feature  is  that  a  normal  simplex- 
type  basis  factorization  B  =  LU  is  suitable  for  the  main  steps  of  the  algorithm.  In 
addition,  large  numbers  of  columns  of  B  can  be  replaced  in  the  basis  at  once. 

8.3.  Shifted  barrier  methods 

Apart  from  Eriksson’s  algorithm,  the  only  method  that  specifically  allows  an  arbi¬ 
trary  starting  point  is  the  shifted  barrier  approach  (Gill  et  al.  [21]).  For  problem 
(2)  the  objective  function  is  of  the  form 

Eu,,,(x)  =  c^i  -  ln(x>  -I-  Sj), 

where  in  is  a  set  of  weights  corresponding  to  /i,  and  s  is  a  set  of  shifts. 

Apart  from  the  ability  to  use  a  “good”  starting  point,  an  advantage  is  that 
the  condition  of  the  system  used  to  obtain  a  search  direction  can  be  controlled 
by  judicious  choice  of  w  and  s.  Convergence  results  have  been  established  and  a 
preliminary  implementation  has  been  developed  at  Stanford. 

9.  Conclusions 

Some  broad  thoughts  have  already  been  given  in  Section  1.4.  Many  additional 
references  could  have  been  cited,  all  contributing  to  the  current  outlook. 
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In  summary,  we  believe  that  the  simplex  method  will  remain  the  workhorse  for 
the  majority  of  existing  applications,  given  practitioners’  normal  mode  of  operation: 
frequent  restarts  on  slightly  modified  models. 

However,  from  the  last  few  years  of  experience  we  would  say  that  interior-point 
algorithms  are  efficient  on  a  rather  larger  class  of  (general)  linear  programs  than 
we  once  would  have  supposed.  New  and  extremely  large  applications  will  require 
decomposition  if  simplex  techniques  are  to  succeed.  The  alternative  is  a  radically 
new  approach,  such  as  we  seem  to  have  here.  While  much  trial  and  error  lies 
ahead,  we  feel  that  the  future  for  the  new  approach  to  linear  programming  shines 
challenging  and  bright. 

In  the  meantime,  we  have  a  new  way  of  identifying  examples  where  the  simplex 
method  has  performed  poorly.  We  also  have  a  renewed  interest  in  improving  the 
simplex  method  (e.g.,  [20,25]). 
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AOMTRACT  (I 

The  world  of  mathematical  programming  has  seen  a  remarkable  surge  of  activity 
following  publication  of  Karmarkar's  projective  algorithm  in  May  1984.  A 
review  of  the  ensuing  three  years  has  already  appeared  (Monma  [31]).  One 
year  later,  we  review  some  of  the  main  methods  and  surrounding  events,  and 
focus  on  references  that  contain  computational  results. 


leiTIOM  OP  I  MOV  M  I*  O' 


MCURltV  CkAMPICATIOM  O 


E  A/  D 

pate 

II' 


